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We present and test a general-purpose code, called PPASPH, for evolving self-gravitating fluids 
in astrophysics, both with and without a collisionless component. In PPASPH, hydrodynamical 
properties are computed by using the SPH (Smoothed Particle Hydrodynamics) method while, 
unlike most previous implementations of SPH, gravitational forces are computed by a PP (Particle- 
Particle) approach. Other important features of this code are: a) PPASPH takes into account the 
contributions of all particles to the gravitational and hydrodynamical forces on any other particle. 
This results in a better energy conservation; b) Smoothing lengths are updated by an iterative 
procedure which ensures an exactly constant number of neighbors around each gas particle, c) 
Cooling processes have been implemented in an integrated form which includes a special treatment 
to avoid a non-physical catastrophic cooling phenomenon. Such a procedure ensures that cooling 
does not limit the timestep. d) Hydrodynamics equations optionally include the correction terms 
(hereafter V/i terms) appearing when h(t,r) is not constant. 

Our code has been implemented by using the data parallel programming model on The Connection 
Machine (CM), which allows for an efficient unification of the SPH and PP methods with costs 
per time step growing as ~ N . 

PPASPH has been applied to study the importance of adaptative smoothing correction terms on 
the entropy conservation. We confirm Hernquist's (1993) interpretation of the entropy violation 
observed in previous SPH simulations as a result of having neglected these terms. An improvement 
on the entropy conservation is not found by merely considering larger numbers of particles or 
different Ns choices. The correct continuum description is only obtained if the V/i correction 
terms are included. Otherwise, the entropy conservation is always rather poor as compared to that 
found for the total energy. 

Subject headings: hydrodynamics - numerical methods 



1. INTRODUCTION 

Most hydrodynamical problems require numerical 
calculations because of their complexity. Several nu- 
merical methods have been developed to solve the 
equations of hydrodynamics, however one of the best 
suited for astrophysical problems is the Smoothed 
Particle Hydrodynamics (SPH) technique. 



SPH is an N-body integration scheme introduced 
by Lucy (1977) and Gingold & Monaghan (1977) as 
an attempt to model continuum physics avoiding the 
limitations of grid-based finite difference methods. 
In SPH, fluid elements constituting the system are 
sampled and represented by particles, and dynami- 
cal equations are obtained from the Lagrangian form 
of the hydrodynamic conservation laws. SPH has two 
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main advantages with respect to other techniques: on 
one hand, since it follows the evolution of individual 
fluid elements, the computational resources are put 
where they are needed most. On the other hand, 
there is no grid constraining the dynamic range or 
the global geometry of the systems being studied. 

Gravitational interactions between particles can 
be obtained from different techniques. The most 
straightforward method consists of computing grav- 
itational accelerations as the direct sum of all in- 
teractions between particles. This approach, called 
particle-particle (PP), has several theoretical advan- 
tages over other potential solvers in the context of 
SPH (Hernquist & Katz 1989). In particular, it is 
fully Lagrangian, it does not use a grid, and energy 
conservation is generally better. However, for simu- 
lations involving a large number of particles, a uni- 
fication of SPH and PP techniques was prohibitively 
expensive because the computing time per step scaled 
as ~ N 2 . Gravitational accelerations in SPH codes 
are then usually computed by using other approaches. 
For instance, grid-based methods (Monaghan & Lat- 
tanzio 1985, Evrard 1988), or the hierarchical tree 
method (Hernquist & Katz 1989). 

Thanks to the development of computers with a 
data parallel programming model as, for example, the 
Connection Machine (CM), we can now envisage the 
use of a PP approach in SPH. On this kind of comput- 
ers, the contribution of a particle to the gravitational 
accelerations of all the other particles can be calcu- 
lated in just one parallel operation (Alimi & Scholl 
1993, Serna, Alimi & Scholl 1994, Scholl & Alimi 
1995). The computing time then scales as ~ N and 
the practical advantages of other potential solvers do 
not hold. 

We will present in this paper a general-purpose 
code, called PPASPH, where SPH and PP techniques 
have been coupled and implemented on CM. After 
exhaustively testing this code, it has been applied to 
analyze in detail the following important aspect of the 
SPH method: 

Although SPH allows for an easy implementation 
of adaptative resolution scales, it introduces addi- 
tional terms in the equations of motion. These ad- 
ditional terms are usually neglected because they are 
computationally expensive and, for a high enough 
number of particles, they are much smaller than 
the other ones (Gingold & Monaghan 1982, Evrard 
1988). However, Hernquist (1993) has found that, 
when adaptative SPH algorithms are used to simu- 



late the evolution of adiabatic systems, a simultane- 
ous good conservation of energy and entropy is not ob- 
tained. If the thermal energy equation is integrated, 
the total entropy is not conserved as accurately as the 
energy. Reciprocally, if an entropy equation is inte- 
grated, then the total energy is not conserved as accu- 
rately as the entropy. Hernquist (1993) then claimed 
that conclusions at high resolution using SPH must 
be accepted with caution. This is in principle also 
extensible to any other fluid algorithm using adapta- 
tive resolution scales or grids. It is then important to 
analyze in detail if this difficulty can in fact be solved 
by including the previously neglected terms or, in the 
opposite, if it reveals some weakness inherent to the 
SPH technique itself. 

This paper is organized as follows. In Section 2 we 
describe our PPASPH code as well as the implemen- 
tation of radiative cooling processes. Several tests on 
this code are then presented in Section 3. The in- 
clusion of additional terms related to the adaptative 
resolution scales is described in Section 4, as well as 
the analysis of their importance on the simultaneous 
conservation of entropy and energy. Section 5 sum- 
marizes our main conclusions. 

2. PPASPH 

2.1. Basic principles of the SPH method 

In SPH, any macroscopic variable (density, pres- 
sure gradient,...), f(r), is conveniently calculated in 
terms of its values at a set of disordered points (the 
particles) by means of an interpolation technique 
known as kernel estimation. This technique is equiv- 
alent to convolving the field f(r) with a smoothing, 
or filter, function W(r — r', h) to produce an estimate 
of the field, fs(f), where local statistical fluctuations 
have been smoothed out: 

fs(r) = J f(r')W(r-r',h)dr' , (1) 

the integration being over all space. The smoothing 
length h specifies the extent of the averaging volume 
and it then determines the local spatial resolution. 
The smoothing kernel W(r — r' , h) is assumed to be 
spherically symmetric, and normalized to unity when 
integrated over the space, as part of the requirement 
lim/^o fs{r) = f{r). 

For numerical work, where a finite number N g of 
gas particles is used, the integral interpolant is ap- 
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proximated by a summation interpolant: 

f s ( ri ) = J2 mj^-Win - rj^ij) , (2) 

j= i py r i> 

where p{r j) is the density at the position of particle 
j. As a particular case of Eq. (2), the smoothed 
estimate of the density at r 8 - is 

P( v i) = m j W ( r ij , hij) . (3) 

Here, r 8 j =| r 8 — Vj |, and hij denotes a sym- 
metrized smoothing length, hij = (hi + hj)/2, nec- 
essary to avoid a violation of the reciprocity principle 
(Evrard 1988). An alternative proposed by Hern- 
quist & Katz (1989) is the use of kernel averages, 
\W(rij,hi) + W(rij,hj)]/2, instead of an average of 
smoothing lengths. Both possibilities provide similar 
accuracies. We have chosen that proposed by Evrard 
(1988) because it requires a slightly smaller amount 
of CPU time. 

The smoothing formalism also provides a natural 
means for estimating gradients of the local fluid prop- 
erties, or any other derivatives. For example, from Eq. 

(1) 

[Vf{r)] s = j Vr>f(r')W(r - r' , h)dr' (4) 

or, after integration by parts and approximation by a 
summation interpolant 

(V/) 8 = £ m^V.ffr,,^) . (5) 

j P( r j) 

2.2. Kernels 

Many kernel functions can be devised in SPH. That 
which allows for an easier physical interpretation of 
any SPH equation is the Gaussian kernel 

where v is the number of dimensions. Although 
this kernel interpolates with high accuracy, it has 
the practical disadvantage that it is not exactly zero 
for finite rij/hij values. Consequently, much parti- 
cles contribute to local properties and, on sequen- 
tial and vectorial computers, much computational ef- 
fort is needed. In order to avoid this inconvenience, 
several authors use instead kernels based on spline 



functions with compact support, as that proposed by 
Monaghan & Lattanzio (1985). 

a { + 1 (o<*o-<i) 

ij { > 2) 

where Sij = f'ijjhij and a is a normalization constant 
with the values 2/3, 10/77T, 1/tt, in one, two, and 
three dimensions, respectively. This kernel has the 
practical advantage that it is exactly vanishing for 

^ij/hij ^ 2. 

In a parallel programming model, the practical ad- 
vantages of spline kernels do not hold because sums 
over all the particles can be performed in a parallel 
way. PPASPH takes then into account all contribu- 
tions to the local properties, even when they are very 
small or exactly vanishing. In all the simulations pre- 
sented in this paper we have used a Gaussian kernel. 
Nevertheless, other possible choices, as Eq. (7), have 
been also implemented in PPASPH. 

2.3. Smoothing Lengths 

Ideally, the individual particle smoothing lengths 
hi must be updated such that each particle inter- 
acts with a constant number of neighbors Ns ■ By 
neighbors we mean those particles j with distances 
r ij ^ Hhi, where % is a constant for each kind 
of kernel 1 . Such a condition can be exactly im- 
plemented without additional computing time when 
gravitational forces are computed using a tree code 
and neighbor lists are then available. In principle, 
other gravitational schemes, as the PP one, would 

— 1/3 

need update hi from conditions like hi oc p i . How- 
ever, this last condition only ensures a roughly con- 
stant number of neighbors and some adjustments are 
needed during a simulation to avoid an excessively 
large deviation from Ns ■ 

We have constructed an efficient algorithm for data 
parallel programming model which exactly updates hi 
while consuming only a modest amount of computing 
time. Since the identity of neighbors is not necessary 
in a PP code, we must just find the sphere centered in 
each particle i which contains a specified number of 
particles. Such an algorithm has been implemented 
in the following way. 

For each SPH particle i, we start with a pre- 
dictor value h" +1 = h" of its smoothing length. 



7i = 2 for a spline kernel, while 7i ~ 3 for a Gaussian kernel 
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We then count the number of particles j such that 
r ij < On CM, this number N t can be fastly 

obtained by means of the CO UNT command. If iVj is 
not equal to Ns , a corrector value is computed from 



n + l 



n + l 1 



(Ns/Ni) 1 ' 3 ] 



(8) 



This corrector value h" +1 is used as predictor h" +1 
for the next iteration step and the sequence is followed 
until that Ni = Ns ■ This procedure can be considered 
as essentially identical to that used by Hernquist & 
Katz (1989), but with a zero tolerance parameter. 

However, in order to avoid iterations with an os- 
cillatory convergence, the predictor hi value given by 
Eq. (8) is forced to be strictly within the interval 
(hj ,h 2 ), where h\ and tif are the last h values found in 
the iteration sequence such that they imply Ni < Ns 
and Ni > Ns , respectively. Since we iterate until 
that Ni = Ns , our method gives by construction an 
exactly fixed number of neighbors for all particles. 
This algorithm is obviously independent of Ns and, 
on a parallel computer, it just grows with the number 
of gas particles as ~ N g . 

The choice of Ns is determined by the condition 
that the theoretical density field at initial conditions 
must be well described by the SPH density field (typ- 
ically, Ns 6 [30, 50]). It must be however noted that, 
if the number of neighbors experiences discrete jumps 
as hi for a particle is increased, our zero tolerance al- 
gorithm could fail for certain choices of Ns ■ Such 
a situation is extremely rare in the course of a sim- 
ulation, but possible at the initial conditions when 
particles are distributed on a regular lattice. In this 
case, the input Ns value must be chosen with cau- 
tion or, otherwise, the number of iterations must be 
limited to a maximum value. 

2.4. Dynamical Equations 

The evolution of particle i is determined by Euler's 
equations 



drj_ 

dt 
dvj 
dt 



Pi 



(9a) 
(9b) 



where $ 8 is the gravitational potential at r 8 , Pi is the 
local pressure and aj lsc is an artificial viscosity term 
allowing for the presence of shock waves in the flow. 



The SPH representation of Eqs. (9) is not unique. 
Several forms of the equations of motion can be found 
in the literature, none of which appears to be clearly 
superior to the others. The most often used SPH ex- 
pression for the pressure gradient and viscosity terms 
is 

V P- 

' + < 1SC = (io) 



Pi 



E ra J 5 + 5 + n 'J V 8 ^(r 8J , %) 
j=\ \ Pi p i J 

where II 8 j stands for the artificial viscosity. 

Several expressions for the artificial viscosity have 
been proposed in the literature. Up to date, that 
giving the best results in SPH codes is (Monaghan & 
Gingold 1983): 



n,- 



Pp'h 



fii) 



Pa 



where a and [3 are constant parameters of order unity, 
Cj is the sound speed at the position of particle i, 
Cij = (ci + cj)/2, Pij = (pi + pj)/2, and 



Pi, 







VijTij < 

VijTij > 



(12) 



with Vij = Vi — Vj and r/ 2 being a softening constant 
avoiding numerical divergences. Typically, r/ 2 = 0.01. 

In order to compute the local pressure, we must 
specify an equation of state Pi(pi,Ui). For an ideal 
gas the equation of state is 



p i = (7 - 1 



(13) 



where 7 = 5/3. Some problems (see section 3.1.2.) 
could however require the isothermal equation of state 
Pi = c? so pi, where q so is the isothermal sound speed. 
In this last case, an equation for the evolution of the 
specific internal energy w 8 - is not needed. 

In general, we used the non-symmetrized thermal 
energy equation: 



dui 
~dt 



j = l ^ P i 



Pi Ih 



Vij'V i W(rij,h i 



•cool 



which gives better results than a symmetrized equa- 
tion for the evolution of w 8 - (Benz 1989). 

The term u co ° l appearing in Eq. (14) denotes 
the cooling (or heating) rates associated with non- 
adiabatic processes other than shocks. 
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2.5. Implementation of Cooling 

The rate of the specific thermal energy variation 
due to radiative cooling processes is 



du 



zool 



dt 



P 



(15) 



where A(T) is the cooling function. This function has 
been obtained by considering the gas as an optically 
thin 'primordial' mixture of H and He (mass fractions 
X = 0.76, Y = 0.24, respectively) in collisional ion- 
ization equilibrium at temperature TJ-. The cooling 
function is then given by A = Ab rem + Ah + Ah 6 , 
where Ab rem is the bremsstrahlung cooling (Tucker 
1975), while Ah and An e are the cooling functions for 
radiative recombination and line emission processes of 
hydrogen and helium (Bond et al. 1984): 

Abrem = 7.31 X 10 22 fT 1 ' 2 



A 



H 



2.32 x 10 25 / 2 



1 + 0.25T 8 



(16) 



A He = 7.01xl0 22 / 2 T3 ( 1 + 3 - 6xl0 " 5T4 ) 



1 + 3.3 x 10" 8 T 8 



/ being the ionization fraction 



/ 



1-/ 



3.2 X 10V13.6/T T 1.22 



(17) 



Temperatures in Eqs. (16)-(17) are written in eV, 
while A is expressed in [erg cm 3 g _2 s -1 ]. 

Usually, the cooling processes impose strong limits 
on the integration timestep. These constraints can 
be however softened by implementing such processes 
in an integrated form controlled by the Courant con- 
dition (see below section 2.8.). As a matter of fact, 
since this last condition ensures that densities do not 
change considerably over one step, equation (15) can 
be integrated to give (Thomas & Couchman 1992): 



du 



zool 



A,- 



At_ 

Pi 



(18) 



Then, we have only to find the Am? 00 ' values which 
satisfy the above equation. 

In some few underresolved zones coincident with 
the recombination front, the previous method could 
overestimate cooling processes. An special treatment 
for these zones must be then included like, for in- 
stance, that proposed by Anninos & Norman (1994). 



Those authors enforce the pressure equilibrium condi- 
tion at the cooling front to avoid such a non-physical 
catastrophic cooling phenomenon. 

We have adapted Anninos & Norman's procedure 
to our SPH code. That is, particles which overcool 
are selected by using two criteria: 



1. The local cooling time (A^ 00 ' = is 
smaller than the viscous-sound crossing time 
At cv (see Eq. [28] bellow): 



At 



cool 



< At c 



(19) 



The local pressure after updating the cooling 
term is smaller than half the average pressure 
Pi of the neighbors of particle i 



Pi < Pi/2 



(20) 



When both criteria are satisfied for a particle, we en- 
force the condition Pi = Pi. 

2.6. Gravitational interactions 

Gravitational interactions in PPASPH are com- 
puted by using the PP method. This approach 
has some theoretical advantages over other potential 
solvers in the context of SPH (see Sect. 1). In addi- 
tion, it is simpler and easier to parallelize on CM than 
other methods as, e.g., the hierarchical tree approach. 

Since particles in SPH just represent elements of 
a continuum fluid, gravitational interactions between 
particles must be smoothed by using the techniques 
of Sect. 2.1. According to the procedure outlined 
by Gingold & Monaghan (1977), for a gravitational 
potential defined as 



-G 



Ptot(r')dr' 



(21) 



the smoothed expression of the gravitational acceler- 
ation of particle i is 



N 



Gm i \ S 



W(r, e)r 2 dr 



(22) 



where a gravitational smoothing length e, different 
from the hydrodynamical one h, must be used. 

The integral appearing in Eqs. (22) can be analyti- 
cally solved provided that the functional form of the 
kernel W has been specified. 
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For the Gaussian kernel given by Eq. (6), we find 



N 



■E 

i=i 

N 

■E 

i=i 



Gnij jerf(sij) 2e-( s ^') 



Gmj erf(sjj) 



r,-i23) 



(24) 



where Sij = rij/cij. The corresponding expressions 
for the spline kernel given in Eq. (7) can be found in 
the Appendix of Hernquist & Katz (1989). 

It must be noted that, when Sij ^> 1, the above ex- 
pressions lead to the Newtonian ones, <f>ij = —Grrij/rij 
and dij = —GrrijVij /rfj while, for s — > 0, they tend 
towards a constant value. For a Gaussian kernel we 
find 4>ij = —2Grrij I eijTT 3 ! 2 , a 8 j = — AGmjTij /'ieijir 1 ! 2 . 
As it has been previously noted by other authors, 
the SPH technique provides a more coherent way 
to softening gravitational interactions that the usual 



2X1/2 



procedure which represents 



■ £y)/2, where the 



particles as Plummer spheres 
In the above equations e 8 -j 
gravitational smoothing lengths do not change with 
time, and constitute a minimum value for the vari- 
able smoothing length h of the gas particles. The 
choice of a minimum value for h is similar to set- 
ting a maximum density value and, indirectly, a mini- 
mum value for the time-step (Navarro & White 1993). 
The justification for such a condition is that it would 
be wasteful to estimate the pressure gradients with 
higher resolution than the gravitational potential in 
regions where the softening is important. 

2.7. Time Stepping 

The time integration in PPASPH is performed 
using a PEC (Predict, Evaluate, Correct) variable 
timestep scheme similar to that considered by Couch- 
man et al. (1994). According to this scheme, one 
enters the time t n with known positions r n , veloci- 
ties v n , and accelerations a' n , for all the N particles, 
as well as the hydrodynamical quantities (smoothing 
lengths h n , specific internal energies u n , and their 
derivatives ii n ) for all the N g gas particles. 

The sequence initiates by predicting variable values 
(denoted by primes) at t n+ i according to 



' n + l 



J n + 1 



-'n + l 



r n +v n At + a' n (At) 2 /2 
v n +a' n At (25) 



The above predicted quantities are then used to 
compute a' n+1 and u' n+ i at r' by using Eqs. (9b) and 
(14). These predicted 'forces' are then used to correct 
the positions, velocities, and thermal energies: 



r„+i 


= r 'n + l ~ 


VA{a' n+1 


-a' n )(At n f/2 


V n + 1 


= v n+l " 




- a' n )At n (26) 


«n + l 


= M n + 1 " 


K^K+i 


- <)A*„ 



where B = 1/2 is required to obtain accurate veloc- 
ities to second order, while the choice of A and C is 
somewhat arbitrary. The choice .4 = gives a scheme 
similar to a leapfrog scheme except that velocities are 
predicted forward to the same time as the positions 
before force evaluation. The choice A = 1/3 nomi- 
nally gives third-order accuracy for positions (but this 
is swamped by the error in velocities). We have cho- 
sen in this paper .4=1/3 and C = 1/2. 

Cooling contributions to the thermal energies are 
finally included by using the integrated scheme of 
Sect. 2.5. 

The advantages of this PEC scheme as compared 
to a Runge-Kutta one have been shown by Couchman 
et al. (1994). They found that, when the timestep 
is slightly overestimated, SPH simulations using the 
PEC scheme remain strongly stable while those using 
a Runge-Kutta one can become exponentially insta- 
ble. 

2.8. Timestep length 

In order to maintain the numerical integration sta- 
bility, the timestep must be modified at each step 
according to different criteria. A first timestep con- 
trol is that concerning the time scale for significant 
displacements or changes in velocity due to accelera- 
tions: 

St a = mm j^J , (27) 

where, for dark-matter particles, the smoothing length 
hi is replaced by the softening parameter e 8 -. 

A second limit on At is usually given by a timestep 
control which combines the Courant and the viscous 
conditions: 



hi 



Ci + 1.2(acj + /? maxj | frj |) 



(28) 



U r . 



it' At 



The time integration of shock contributions to 
the thermal energy equation is also limited by the 
Courant condition. However, if the cooling function 
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A is nonzero, an additional timestep control involv- 
ing m coo i should be needed and would imply very se- 
vere limits on At. Since cooling processes in PPASPH 
are included in an integrated form with the enforced 
pressure equilibrium treatment described in Sect. 2.5, 
such a timestep control is not necessary. 
The timestep is then given by 

At = mm(8t a ,8t cv ) . (29) 

If the At resulting from Eq. (29) is greater than 
some input value (At) max , the timestep is forced to 
be equal to this upper limit. 

3. TESTS OF PPASPH 

PPASPH has been applied to a number of systems 
in order to test its ability to reproduce known analytic 
or numerical solutions. 

3.1. One-dimensional tests 

All ID-tests presented here were performed with 
N = 4096, N s = 40, a = /? = 1 and rj 2 = 0.01. Dissi- 
pational effects, other than those associated with the 
artificial viscosity, were ignored, as well as gravita- 
tional interactions. 

3.1.1. The shock tube problem 

The one-dimensional shock tube problem proposed 
by Sod (1978) has become a standard test of all trans- 
port and source terms (including artificial viscosities) 
of hydrodynamic algorithms. It considers a perfect 
gas distributed on the «-axis. A diaphragm at xo ini- 
tially separates two regions which have different den- 
sities and pressures. All particles are initially at rest. 
At time t = the diaphragm is broken and both re- 
gions start to interact. Nonlinear waves are then gen- 
erated at the discontinuity and propagate into each 
region: a shock wave which moves from the high to 
the small pressure region, while the associated rar- 
efaction wave moves in the inverse sense. At the 
contact discontinuity, the fluid density and specific 
energy are discontinuous, while the velocity and pres- 
sure are continuous. However, in the location of the 
shock wave, all quantities (P, p, v and u) will be dis- 
continuous. The analytical solution to this problem 
has been given by Hawley et al. (1984), and Rasio & 
Shapiro (1991). 

In our simulation, we have initially considered a 
7 = 1.4 perfect gas distributed in the interval — 1 < 



x < 1 according to: 

p=lP=lv = (for x < 0) 

p = 0.25 P = 0.1795 v = (for x > 0) 

Fig. 1 shows our results at t=0.15. We see from 
this figure that our results are in excellent agreement 
with the analytical solutions. The resulting profiles 
both in the shock wave (located between x ~ 0.2 and 
x ~ 0.25) and in the contact discontinuity (located at 
x ~ 0.1) are much less rounded than in previous SPH 
computations (see, e.g., Monaghan & Gingold 1983, 
Hernquist & Katz 1989, Rasio & Shapiro 1991) as a 
result of having used a larger number of particles and, 
hence, a better resolution. 

Much more encouraging for us is the almost com- 
plete suppression of postshock oscillations in our re- 
sults. These oscillations can be seen in the previ- 
ous SPH simulations of this problem, especially in 
the velocity field, while no high-frequency vibrations 
are perceptible in our results. Moreover, the over- 
shoot observed in the velocity field at the tail of the 
rarefaction wave (located at x ~ —0.05) is smaller in 
our results than in those previously obtained in the 
literature. 

The weak blip observed in the pressure profile at 
the contact discontinuity is normal in SPH codes. 
Such non-physical blip has been explained by Mon- 
aghan and Gingold (1983) as due to the fact that the 
smoothed estimate of pressure is computed by using 
discontinuous quantities. It is then inevitable that P 
has some slight perturbation at the contact disconti- 
nuity, but it has a negligible effect on the motion. 

3.1.2. The isothermal shock problem 

Another ID-problem often used to test hydrody- 
namics codes is that proposed by Leboeuf, Tajima & 
Dawson (1979). This problem initially considers an 
isothermal fluid with a square wave type density pro- 
file consisting of a central dense region of density ph , 
surrounded by dilute regions with density pi < ph- 
The fluid evolution is characterized by the rapid for- 
mation of a density plateau between rarefaction and 
shock front regions. The shock front speed, v s , as 
well as the density p p and the velocity fluid Vf in the 
plateau region can be analytically obtained as a func- 
tion of Ph/ Pi by using the solutions given by Leboeuf, 
Tajima & Dawson (1979). 

In our SPH simulation we have considered a initial 
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Fig. 1. — a) Density, b) pressure, c)velocity, and d) specific internal energy profiles at t=0. 15 in the one-dimensional 
shock tube problem. Points represent the PPASPH results, and solid lines are the analytical solutions. 
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density profile given by 



x < 7 or x > 7 p = 4 
-7 < x < 7 p = 9 

where units were chosen so that c 2 = 1. 

This profile implies the theoretical values p p = Q 
and = 0.4. 

The results given by PPASPH are very satisfactory 
as compared to those previously obtained by other au- 
thors (e.g., Gingold & Monaghan 1982, and Haddad, 
Clausset & Combes 1991) using the SPH method. 
The theoretical values for p p and Vf are accurately 
obtained in our simulations (Fig. 2). Moreover, the 
plateau regions are almost free from high frequency 
postshock oscillations and broadening effects, which 
were instead very present in the results reported by 
the above quoted authors. An overshoot is however 
observed in the density profile at each corner of the 
high density region. No other similar overshoots are 
found in the other discontinuities both in the density 
and in the velocity profiles. 




3.2. Adiabatic collapse of a non-rotating gas 
sphere 

A 3D-problem usually considered to test hydrody- 
namical codes is that concerning the adiabatic col- 
lapse of a non-rotating gas sphere. This problem 
has been studied from a finite-difference method by 
Thomas (1987), and from SPH simulations by Evrard 
(1988) and Hernquist & Katz (1989). In order to facil- 
itate the comparison of our results to those obtained 
by these authors, we have taken their same initial con- 
ditions, that is, we initially consider a gas sphere of 
radius R and total mass Mt , with density profile 



P 



2irR 2 



1 



(30) 



All the N = 4096 gas particles have initially the same 
specific internal energy u = O.O&GMt/ R. The ratio of 
specific heats is 7 = 5/3 and the Gingold & Monaghan 
artificial viscosity (Eq. 11) was used with a = 1, 
(3 = 2 and r/ 2 = 0.01. The gravitational softening 
parameter was e = 0.05 and a Gaussian kernel was 
assumed with Ns = 40. Units were taken so that 
G = M T = R= 1. 

The time evolution of different system profiles are 
shown in Figs. 3 at times matching those given by 
Evrard (1988) and Hernquist & Katz (1989). Our re- 
sults are in excellent agreement with those presented 




Fig. 2. — a) Density, and b) velocity profiles at t=4 in 
the one-dimensional isothermal shock problem. Units 
are defined by c 2 = 1 
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by those authors. Initially far from equilibrium, the 
system collapses converting most of its kinetic energy 
into heat (between t « 0.8 and t « 1.3). A slow ex- 
pansion follows and, at late times, a core-halo struc- 
ture develops with nearly isothermal inner regions and 
the outer regions cooling adiabatically. The location 
and strength of the shock are also well reproduced by 
our code. For example, we find that, at t = 0.88, the 
shock appears located at r « 0.2 with a Mach num- 
ber of « 4. Viscosity erases very efficiently radial mo- 
tions in the central regions while the collapsing outer 
regions are rebound and, later (t ;> 1.3), they expand 
constituting a rarefied, adiabatically cooling halo. In- 
ternal regions at t = 0.88 have M-profiles slightly in- 
creasing with r but, at the end of our simulation, these 
regions have a nearly constant temperature. This last 
feature agrees better with Thomas's finite-difference 
results than the previously quoted SPH simulations, 
where a small rise in the thermal energy profile was 
still present in the final configuration. 

Other particular features also agree with those 
obtained by Evrard (1988) and Hernquist & Katz 
(1989). For example, the relative maximum observed 
at t = 2.2 in the v/c profile which, although less ev- 
ident in our simulation, is located at r « 0.6. Obvi- 
ously, there exist very small differences between the 
results reported by all these codes. They just come 
from the different resolutions and artificial viscosity 
expressions used in these simulations. 

3.3. Test of Cooling Simulations 

In order to test simulations including radiative 
cooling processes, we have simulated the collapse of 
a rotating sphere ('protogalaxy') with a dominant 
amount of dark matter. This numerical experiment 
can be compared to that performed by Navarro & 
White (1993). 

Initial conditions consist of N gas = Ndm = 1736 
particles. Positions are obtained from a p(r) a r" 1 
spherical perturbed grid, while velocities are chosen so 
that the sphere will be in solid-body rotation around 
the z-axis with a spin parameter A = J \ E I 1 / 2 
/GM ( 5 /( 2 ~ 0.1 (J and E stand for the total angu- 
lar momentum and total energy, respectively). The 
initial radius of the sphere was R tot = 100 kpc, and 
its total mass was M to t = 10 12 Mq. The gas repre- 
sents a 10 per cent of this mass, and is initially at a 
uniform temperature, T = 10 3 K. Gravitational soft- 
ening parameters were taken to be 2 and 5 kpc for 



the gas and dark matter, respectively, and units were 
chosen so that G=1,[M] = 1O 1O M , [L] = 1 kpc. 

When radiative cooling processes are switched-on, 
the thermal energy gained by particles through shocks 
is quickly radiated away. Consequently, the gas never 
gets heated to the virial temperature (ss 3 x 10 6 K). 
Without pressure support, collapse of the gas pro- 
ceeds unimpeded until it becomes centrifugally sup- 
ported in a thin disc-like structure (Fall & Efstathiou 
1980, White 1991). As in Navarro & White's (1993) 
simulation, we find that the disc is almost completely 
formed few after the collapse time (t ~ 120) and 
evolves little thereafter (see Figs. 4). Shocks dissi- 
pate very efficiently the energy in radial motions and 
the remaining kinetic energy is invested in the rota- 
tional motions that support the disc. 

A spiral-like structure in the gaseous disc starts to 
be apparent at t « 160, and remains during the rest of 
the simulation. The aspect of this spiral structure at 
t = 320 is intermediate between the two simulations 
reported by Navarro & White (1993). These authors 
considered two cases, with gas mass fractions of 0.1 
and 0.02, respectively. In the first case, they found a 
locally unstable disc broken into small clumps while, 
in the second case, they found a much cleaner spiral 
structure. Although we used the same gas mass frac- 
tion than in the first of such experiences, we find a 
more locally stable disc. This is certainly due to dif- 
ferences in the initial conditions. The sound speed or 
the circular frequency are probably higher in our sim- 
ulation and, consequently, Toomre's (1964) stability 
parameter is also higher than in the first simulation 
by Navarro & White (1993). Other features are in- 
stead similar to those obtained by the above quoted 
authors. For example, the presence at t = 320 of two 
or three main spiral arms, and a dense and small core 
surrounded by a more dilute region. 

4. CORRECTIONS FOR ADAPTATIVE - 
SMOOTHING LENGTHS 

We have previously neglected in Eqs. (9b) and 
(14) the terms resulting from the fact that h is not 
a constant. As quoted in the introduction, for not 
very high numbers of particles, such approximation 
could introduce non-negligible errors on the SPH re- 
sults. Hernquist (1993) claimed that the poor entropy 
conservation observed in some SPH simulations is due 
to having neglected such terms. Although the global 
properties of a system do not seem to be altered, ex- 
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Fig. 3. — Density, pressure, specific internal energy, radial velocity, and Mach number profiles in the adiabatic 
collapse of a non-rotating gas sphere. The times shown are indicated in the upper right corner of each frame. Units 
are G = M T = R = 1. 
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cept for a very slight delay in the collapse time of 
structures, conclusions involving high resolution as- 
pects should be accepted with caution. 

In this section, we will describe how the terms 
resulting from an adaptative smoothing, called V/i 
terms, have been included in PPASPH, and we will 
analyze their influence on the entropy conservation. 

4.1. The V/i correction terms 

The general expressions for the V/i correction 
terms have been given by Nelson & Papaloizou (1993) 
in the case of smoothing lengths symmetrized as 
(hi + hj)/2 and computed from a procedure of the 
same type than that described in Sect. 2.3. After 
some straightforward algebra, Nelson & Papaloizou's 
expressions can be written in a more compact form 
as: 




m p] dhij [ 



where a 8 j and the correction terms to be added 

in Eqs. (9b) and (14), respectively. Subscripts k m 
denote the most distant neighbor of particle k, that 
is, that particle satisfying \ rk — rk m |= Hh^- 

Notice that the above equations require that the 
most distant neighbor of each particle k be identi- 
fied. We have implemented it by using the CM func- 
tion MAXLOC(rjfe , mask) which locates the maxi- 
mum element, k m , of the array containing the dis- 
tances to k of all particles satisfying the condition 
mask = rjk < Hh^- Although this is performed in a 
parallel way, the function MAXLOC involves commu- 
nication between processors and, therefore, it breaks 
somewhat the parallel efficiency of PPASPH. Com- 
puting times per step in simulations including V/i 
corrections are typically longer by a factor of two than 
those neglecting such terms. 

4.2. Influence of V/i correction terms 

Simulations as that shown in Sect. 3.2 lead at late 
times, t ;> 3, to an equilibrium sphere with the density 
profile displayed in Fig. 3. This model is close enough 
to the collapsed systems which should be found in 



simulations studying the formation of structures. In 
order to analyze the influence of the V/i correction 
terms on the entropy conservation, we have simulated 
the adiabatic evolution of such kind of spheres. 

Initial conditions were then generated by perform- 
ing a simulation like that described in Sect. 3.2, 
but for different numbers of particles. At t = 3, we 
switched-off its self-gravity and viscous pressures (by 
setting a = [3 = 0) in order to ensure that the subse- 
quent evolution must conserve the total entropy. In 
absence of gravitational interactions, this system ex- 
pands fastly and, at t = 3.3, its central density has de- 
creased by a factor of « 25. The evolution from t = 3 
tot = 3.3 must conserve both the total energy, E, and 
the total entropy variable, St = m j \°g[ a j( s )], 
where the entropic function a(s) is defined by 



7 



a{s)i = 

Pi Pi 



1 



(33) 



We have run a series of simulations, where all models 
started from the same type of initial conditions. Their 
evolution was always followed by integrating the en- 
ergy equation (14). The relative variation of energy 
and entropy from t = 3tot = 3.3is shown in Table 1 
(runs labeled by A). We see that, when the V/i cor- 
rection terms are neglected, energy is conserved very 
accurately but there exists a considerable violation in 
the total entropy variable (about 5% in the considered 
time interval). In the opposite, when such correction 
terms are taken into account, both total energy and 
total entropy are conserved very accurately (about 
0.02%). We thus confirm Hernquist's (1993) interpre- 
tation of this entropy violation as a result of having 
neglected these correction terms. 

Inspection of Table 1 shows moreover that the im- 
portance of having neglected the V/i correction terms 
does not decrease when a larger number of particles 
is considered. The entropy violation in these simula- 
tions is in fact about 5% whatever the value of N is. 
However, it must be noted that just taking N — > oo 
and /z — )> is not enough to obtain the proper con- 
tinuum limit. The condition Ns — > oo is also nec- 
essary to ensure that the discrete SPH approach be- 
comes close enough to the continuum physics. Since 
h (x (Ns/N) 1 ' 3 , the joint limit N oo, h and 
Ns —7- oo can be reached by taking Ns oc N a with 
< a < 1. 

In order to analyze more in detail how the lack of 
entropy conservation depends on N and Ns, we have 
performed some further simulations where Ns was 
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taken as Ns oc y N . In some first test experiences, 
with initial conditions generated as before, the AS 
results for different choices of Ns appeared rather in- 
fluenced by differences in the initial conditions and in 
the initial small-scale random noise. Error bars in the 
entropy violation were then rather large (AS ~ 5±2% 
for the considered time interval) and, hence, it was 
difficult to extract reliable conclusions. A possible Ns 
dependence could be masked within the error bars. In 
order to reduce these initial random fluctuations, we 
have considered the outcome of the simulation pre- 
sented in Fig. 3 as that giving the theoretical profiles 
of a collapsed non-rotating sphere. Particle positions 
were then settled according to the theoretical density 
profile by means of the procedure proposed by Whit- 
worth et al. (1995). 

The energy and entropy violation found in these 
numerical experiences is shown in Table 1 (simula- 
tions labeled by B). We see again that, if the V/i 
correction terms are excluded, the outcoming entropy 
violation is much larger than that found for the to- 
tal energy. The lack of entropy conservation does not 
exhibit any systematic dependence on N or Ns ■ Con- 
sequently, if these correction terms are neglected, the 
SPH results never converge towards those implied by 
the correct continuum physics, at least in that con- 
cerning the entropy. 



TABLE 1 

Entropy and Energy Conservation for 
different 



ijiven by 



Run 


V/i terms 


N 


N s 


AE 


AS 


Al 


Excluded 


1024 


40 


0.01% 


4.8% 


A2 


Included 


1024 


40 


0.02% 


0.02% 


A3 


Excluded 


2048 


40 


0.01% 


5.1% 


A4 


Included 


2048 


40 


0.02% 


0.01% 


A5 


Excluded 


4096 


40 


0.02% 


5.3% 


A6 


Included 


4096 


40 


0.02% 


0.02% 


Bl 


Excluded 


2048 


45 


0.02% 


4.7% 


B2 


Included 


2048 


45 


0.03% 


0.02% 


B3 


Excluded 


4096 


64 


0.02% 


5.0% 


B4 


Included 


4096 


64 


0.02% 


0.03% 


B5 


Excluded 


8192 


90 


0.01% 


4.6% 


B6 


Included 


8192 


90 


0.02% 


0.02% 



dt 



j 



(34) 



where p i denotes the sum of adaptative smoothing 
corrections terms on dpi/dt. Differently from the 
terms u 8 and a 8 , the p i term cannot be switched- 
off when a SPH code as that described in Sect. 2 
is used. As a matter of fact, since densities are esti- 
mated in practice by evaluating Eq. (3) at each time 
step, rather than by integrating dpi/dt, the p i term 
is always implicitly incorporated. 

By differentiating a(s) 8 (Eq. [33]) with respect to t, 
we find that the entropic function of particle i changes 
as 

1 da,i _ , _ ^ ft 

di dt Ui pi 

where we have replaced dui/dt and i 
expressions (14) and (34). 

According to Hernquist's (1993) interpretation, if 
Ui is switched-off, the right hand side of Eq. (35) 
only contains the term in p i and, consequently, any 
variation of pi should imply a 'non-physical' variation 
of entropy. In the opposite, if u 8 is switched-on, it 
balances the contribution of p i and, in the absence 
of dissipation, a(s) should be constant particle by 
particle. The results of our simulations show in fact 
that both terms balance very accurately. Therefore, 
Uijui = (7 — l)Pi/pi and, if the adaptative correction 
terms are switched-off, Eq. (35) becomes 



(35) 

i/dt by their 



1 dcii 1 ^-^ 

dt ~ 11: ^ ' 



with 



dhjj 
dr-i 



1 

m 



Pi dWjj dh^ 

' p'l dhij dr-i 



' JJr 



(36) 



(37) 



In order to estimate the importance of terms ap- 
pearing on the right hand side of Eq. (36), we can 
follow the same kind of simple arguments than those 
considered by Evrard (1988). That is, we approxi- 
mate dhij/dri ~ hij /% and, on the other side, we 
take into account that for a Gaussian kernel: 



dWjj 

dhij 



_3_ 

hij 



3 \h{j 



(38) 



Our numerical results can be interpreted as follows: 
If pi is computed from Eq. (3), its time variation is 



where, since the dominant contribution to the hydro- 
dynamics of any particle comes from scales r 8 j « hij , 
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we can expect that the average of [3 — 2(rij/hij) 2 ] is of 
order unity. Equation (36) can be then approximated 
to 

a% dt ~ nu iP ^ mjWij - Pi ' 

and the total entropy variation is nearly given by 

dt ^ Uui pi y ' 

i 

In this equation, the only quantity which depends on 
the number of particles is m 8 - oc N -1 . The contribu- 
tion of each particle to dS/dt then decreases as TV -1 . 
However, since the sum is performed over the N par- 
ticles, the total entropy variation does not depend on 
the number of gas particles. 

5. SUMMARY 

A general-purpose code (PPASPH) for evolving 
self-gravitating fluids in astrophysics, both with and 
without a collisionless component has been described. 

In PPASPH hydrodynamical properties are com- 
puted by using the SPH (Smoothed Particle Hydrody- 
namics) method, while gravitational forces are com- 
puted by a PP (Particle-Particle) approach. A uni- 
fication of the SPH and PP techniques has several 
advantages. Since both techniques are gridless, the 
resulting code is fully Lagrangian and without limita- 
tions on the system geometry, or mesh-related limita- 
tions on the dynamic range in spatial resolution. The 
energy conservation is also generally better. 

This code has been implemented on the massively 
parallel computer The Connection Machine (CM), 
which allows for an efficient unification of the SPH 
and PP methods with costs per time step growing as 
~ N. Moreover, on CM it is also possible to take 
into account, with a minimal cost in computing time, 
the contribution of all particles to the local properties 
of any other particle. Smoothing lengths can also 
be updated so that they imply an exactly constant 
number of neighbors around each particle. 

PPASPH has been applied in order to study the 
importance of correction terms related to adaptative 
algorithms. We have found that the poor entropy 
conservation observed in adaptative SPH simulations 
can in fact be completely improved by taking into ac- 
count the V/i correction terms. In that case, both 
energy and entropy are conserved with the same de- 
gree of accuracy. We thus confirm Hernquist's (1993) 



interpretation of the entropy violation as a result of 
having neglected such correction terms. An improve- 
ment on the entropy conservation cannot be found 
by just taking a larger number of particles and/or a 
larger Ns value. The correct continuum description 
is only obtained when the V/i correction terms are 
included. Otherwise, the entropy conservation is al- 
ways rather poor as compared to that found for the 
total energy. 

We thus conclude that SPH conclusions concern- 
ing high resolution aspects must be accepted with 
caution, even if a very high number of particles has 
been used. For this kind of analyses we believe neces- 
sary to perform some additional simulations including 
the V/i correction terms in order to verify that con- 
clusions are not altered. In principle, these cautions 
ought to apply also to adaptative grid codes. The fact 
that one can evaluate the correction terms explicitly 
in SPH constitutes an additional advantage of this 
approach over grid-based methods. 
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(CEA), Bruyeres le Chatel, France. N-body compu- 
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National de Calcul Parallele en Sciences de la Terre, 
Paris, France. A.S. thanks also the Ministerio de Ed- 
ucacion y Ciencia, Spain, for a post-doctoral fellow- 
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APPENDIX A. 

IMPLEMENTATION ON CM AND TIMING 
ANALYSIS 

In a simulation using particles, the computation- 
ally most expensive parts are those containing a dou- 
ble loop over particles. This is the case in SPH of any 
smoothed estimate and of gravitational accelerations. 

On computers with a data parallel programming 
model, as the Connection Machine, one physical or 
virtual processor is assigned to each particle. The 
contribution of each particle to the smoothed estimate 
/ of all the other particles can be performed by fol- 
lowing a direct summation approach with a comput- 
ing time which grows as N . In such approach, mf/p 
is stored as a sequential array F while kernels and dis- 
tances from particle j to all others are computed and 
stored as a parallel array W(:). The contribution of j 
to the smoothed estimate / of all the other particles is 
then fj(:) = F(j)*W(:). A single loop over j calculat- 
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ing /(:) = /(:) + fj (:) will then lead to the / values for 
all particles. The storage of some quantities both in 
parallel and in sequential arrays, although not strictly 
necessary, results in substantial speed gains by avoid- 
ing communication between processors. Such storage 
can be performed very efficiently by using specific CM 
subroutines. 

The above algorithm then computes all interac- 
tions disregarding the fact that much of them are al- 
most or exactly vanishing. Other algorithms could 
be constructed as, for instance, by do not comput- 
ing contributions when rij/hij > %, or by perform- 
ing the loop over j only for the N s neighbors of each 
particle. However, since CM is conceived to perform 
dummy parallel operations, the first of the above al- 
ternatives would need nearly the same time than our 
algorithm and, in that concerning the second possibil- 
ity, it would be less efficient on CM because it requires 
the individual identification of the neighbors of each 
particle. 

TABLE 2 

CPU times in PPASPH (on CM-5 32 proc.) 



Section of Code CPU time Percentage 



Updating h 


1.85 


5.6 


Hydrodynamics 


4.32 


13.1 


Gravitation 


12.06 


36.7 


Radiative Cooling 


0.56 


1.7 


V/i terms 


14.09 


42.8 


Miscellaneous 


0.03 


0.1 



For illustrative purposes, we show in Table 2 an ex- 
ample of the CPU distribution in the current version 
of PPASPH. This example corresponds to the case 
of the collapse of a N = 4096 non-rotating sphere 
like that considered in Section 3.2, but with radiative 
cooling processes and V/i correction terms switched- 
on. 
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